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A paradigm based on the absolute equilibrium of Galerkin-truncated inviscid systems to aid in understanding turbulence 
[T.-D. Lee, "On some statistical properties of hydrodynamical and magnetohydrodynamical fields," Q. Appl. Math. 10, 
69 (1952)] is taken to study gyrokinetic plasma turbulence: A finite set of Fourier modes of the collisionless gyrokinetic 
equations are kept and the statistical equilibria are calculated; possible implications for plasma turbulence in various 
situations are discussed. For the case of two spatial and one velocity dimension, in the calculation with discretization 
also of velocity v with N grid points (where N + 1 quantities are conserved, corresponding to an energy invariant and 
N entropy-related invariants), the negative temperature states, corresponding to the condensation of the generalized 
energy into the lowest modes, are found. This indicates a generic feature of inverse energy cascade. Comparisons are 
made with some classical results, such as those of Charney-Hasegawa-Mima in the cold-ion limit. There is a universal 
shape for statistical equilibrium of gyrokinetics in three spatial and two velocity dimensions with just one conserved 
quantity. Possible physical relevance to turbulence, such as ITG zonal flows, and to a critical balance hypothesis are 
also discussed. 



[. INTRODUCTION 

Plasma dynamics encompasses a hierarchy of scales with 
distinct physical processes. At scales much larger than the 
mean free path and gyroradius, and time scales much larger 
than the collision time and gyroperiod, the magnetohydrody- 
namics (MHD) model is good (and is often quite useful over 
a wider range of collisionality, particularly for phenomena 
where the parallel kinetic dynamics are not important); while, 
in the opposite limit of high frequencies and small scales, a 
complete kinetic description with the Boltzmann or Vlasov 
equation is necessary. In between, for frequencies well below 
the ion cyclotron frequency but that may still involve scales 
comparable to the gyroradius, a detail of the particle helical 
motion around the field line, the cyclotron angle, may be aver- 
aged out, resulting in a reduced system called gyrokinetics 
With one dimension (the cyclotron angle) and the fast time 
scales associated with that dimension excluded, gyrokinetics 
helps the tractability of turbulent kinetic cascades of plasma 
turbulence numerically and analytically. 

In this contribution, we will present the equilibrium statis- 
tical mechanics of the Fourier Galerkin truncated gyrokinetic 
system and discuss the possible implications for plasma tur- 
bulence. 

Equilibrium-statistical-mechanics approaches to explore 
turbulence have long been attempted to identify the flows or 
to provide some relevant solutions to track the mechanisms 
of fluid turbulent motions, which have been very illuminat- 
ing and promising, if not completely successful!^ One sim- 
ple but efficient strategy, initiated by Leej^ is calculating the 
Gibbs statistics of the Galerkin-truncated system: The flow 
of the Euler equation in phase space is incompressible (where 
the coordinate axes <7j (k) of this phase space are the real and 
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imaginary parts of the Fourier amplitude of the incompress- 
ible velocity field with an upper bound of the wave number 
k), i.e., the dynamics of tJi(k) satisfies the Liouville theo- 
rem, by which an equipartition of energy, which was consid- 
ered as the conserved quantity, among as was then predicted 
(c.f. Appendix |A| for a pedagogical elaboration). There are 
several reasons that the study of the statistical mechanics of 
such idealized systems can be of interest.E2M2l First is that this 
can give analytic (or semi-analytic) predictions for the equi- 
librium statistics that can be used as a nonlinear benchmark 
to test codes. Such nonlinear analytic tests are rare and thus 
valuable. (This has been useful for fluid codes and, in plasma 
physics, for particle-in-cell codesfSKl and could be used for 
continuum kinetic codes as well.) Second, such analytic spec- 
tra can also be useful test cases for analytic theories of tur- 
bulence. Equilibrium statistics has been shown to have subtle 
and deep relevance to statistically nonequilibrium turbulence. 
It has been used to provide insights into two-dimension al (2D) 
guiding-center plasma and 2 -D vortex fluid models,^ ! 18 ! 19 !, 
and other plasma modelPEU. More recently, it has provided 
insights 22 into the unexpected phenomena of spontan eous 
"spin-up" in bounded 2-D fluid turbulence simulations! 2 ^ 2 ^ 
(Interestingly, a current research topic in the fusion field is 
spontaneous rotation observed in tokamaks. 25 26 ) The most 
well-known result from this approach may be the prediction 
of inverse energy cascade in two dimensional turbulence by 
Kraichnanp2l following which Frisch et alP^ calculated the 
magnetohydrodynamic (MHD) absolute equilibrium and il- 
lustrated how the inverse cascade of magnetic helicity may 
help explain the generation of large-scale magnetic fields in 
some astrophysical systems. Another example is how the con- 
cept of 'partial thermalization' has recently been used to un- 
derstand some observed phenomena such as the 'bottleneck' 
near dissipation scales in Fourier space and the reduction of 
intermittency, or its scaling, in physical space, 29 30 which em- 
phasizes the persistence of some aspects of equilibrium statis- 
tical mechanics in turbulence, complementing the other side 
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of our knowledge of the persistence of aspects of cascade 
physics beyond the inertial range (see, e.g., ZhvF^. Revis- 
iting and further extending such powerful tools to accumulate 
relevant knowledge and to examine the relevance to definite 
realities is then important. More recently, this approach has 
been taken to analyze Hall MHD by Servidio et al.J^l find- 
ing that, among others, equipartition of kinetic and magnetic 
energy predicted by Lee^ for Alfvenic MHD turbulence no 
longer holds. Here we will take this paradigm to investigate 
the gyrokinetic model of plasma turbulence. The nontrivial 
new feature in our problem is that the integrations over the 
distributions are functional integrals because of the extra de- 
pendence on velocity of the gyrokinetic variable. 

More generally, understanding the statistical mechanics of 
truncated gyrokinetics can help shed light onto the general na- 
ture of nonlinear coupling in these equations, and phenom- 
ena such as direct or inverse cascades. A better understand- 
ing of nonlinear processes in gyrokinetics may also help in 
the development of more effective sub-grid models for Large- 
Eddy Simulations, and could improve understanding of the ul- 
timate heating mechanisms as the fluctuations cascade to very 
small spacial and velocity scales where collisional dissipation 
occurspEH 

Related to these sub-grid dissipation issues, the three- 
dimensional (3D) energy spectrum of thermal fluctuations that 
we calculate here for a discretized Eulerian gyrokinetic al- 
gorithm turns out to be closely related to the noise spectrum 
calculated earlier for particle-in-cell (Lagrangian) gyrokinetic 
algorithms.^ 

In the two-spatial-and-one-velocity-dimension case, the 
negative temperature state, leading to the condensation of the 
generalized energy at the lowest modes, indicates a generic 
feature of inverse energy cascade. Comparisons are made with 
some classical results, such as those of Charney-Hasegawa- 
Mima in the cold-ion limit, though more generally the spectra 
are modified by finite Larmor radius (FLR) effects which de- 
pend on the temperature parameters. The shape of the statis- 
tical equilibrium for gyrokinetics in three spatial and two ve- 
locity dimensions, where there is just one conserved quantity, 
has a universal energy spectrum shape, resulting from FLR 
effects. 

In the main body we emphasize the general conceptual 
ideas with only necessary details for illustration; specific 
mathematical calculations, physical examples and other inter- 
esting digressions are referred to the appendixes for further 
interests. 



II. FORMULATING THE PROBLEM AND CALCULATING 
THE ABSOLUTE EQUILIBRIA 

To be self-contained, here we very briefly introduce the 
nonlinear gyrokinetic theoretical framework under which we 
will be working. We won't review the complete history of 
the linear and nonlinear gyrokinetic theories^ but will just 
present the basic ideas and results, borrowing from some of 
the treatment and notation of Plunk et alP^U There are sev- 
eral published derivations of gyrokinetic equations with varied 



assumptions and techniques, including recent papers with a 
tutorial emphasis. 5 6 The starting point is the Boltzmann equa- 
tion for the particle distribution function f a (r, v, t) for plasma 
species s located at r moving with velocity v at time t: 
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Here the operator C[f] accounts for the effects of collisions 
and the particles with mass m s and charge q s are accelerated 
by the electric (E) and magnetic (B) fields, which are subject 
to the classical Maxwell equations. The next step is introduc- 
ing the gyrokinetic ordering (which is fundamentally to focus 
on fluctuations that are low frequency compared to the fast 
gyromotion of particles around the magnetic field) and the re- 
sulting expansion parameter. A key operation in the resulting 
equations is the average of any particular quantity around 
a ring of gyroradius p perpendicular (_L) to the magnetic field 
direction (||) surrounding the gyrocenter rP 
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Using a Fourier representation ^(r) = J^k exp(— ik • r)^, 
and considering a straight magnetic field for simplicity here, 

this becomes (^)r = J\ ex P( — ^ ' R)Jo{k±p)'&k, where Jo 
is a Bessel function. 

Writing v = vj_ + vni and / = Fq\ + h + h.o.t. (and 
suppressing the species subscript s for now), with h.o.t. rep- 
resenting "higher order terms," the resulting gyrokinetic equa- 
tions for the case of slab geometry with a homogeneous 
plasma in a straight equilibrium magnetic field Bo = Bqz) 
is 
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complemented with the similar ordering-gyroaveraging treat- 
ment of the Maxwell equations for the electrostatic potential 
ip and the perturbed vector potential A which compose the gy- 
rokinetic potential \ = tp — v • A/c. Here the collisional term 
is omitted. The zeroth and first order term Foi is in general 
taken to be the equilibrium Maxwell distribution (Fq) multi- 
plied by a Boltzmann factor, exp(— g< g/7n ) w 1 — qip/To. In 
what follows below, as in Plunk et al.j 35 ! 36 ! we will work with 
the gyroaveraged, perturbed, guiding center distribution func- 
tion g = h — Foq(cp) n/To, instead of with the non-adiabatic 
component h, and for simplicity we will focus on the case 
of electrostatic fluctuations (neglecting magnetic fluctuations, 
A = 0) with one particle species governed by the gyrokinetic 
equation and the other species having a Boltzmann response 
of some form (discussed below). 

To make it easier to compare with other codes and theo- 
ries that use a variety of normalizations, and in particular to 
make it easier to take the cold-ion limit in 2-D to compare 
with the Hasegawa-Mima equations, we will use a generalized 
normalization for space and time scales based on a reference 
temperature T r , a reference sound speed c r = ^T r /m, and a 
reference gyroradius p r = c r /fi c , (here the mass m and Lar- 
mor (cyclotron) frequency 51 c = qB/mc are for the species 
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that is governed by the gyrokinetic equation), but still scale v {l 
and the velocity dependence of F and g to = WTo/m, 
where To is the temperature of the gyrokinetic species. More 
specifically, we use the following normalizations and defini- 
tions, with physical (dimensional) variables having subscript 
P": 

t = t p c r /L x = Xp/p r y = y P /p r Z = Z p /L 

V - V±AUV lO-lO h-b V&iL Fn - F0pV±3 

X 'U ~~ V th V - T rPr 11 - n P n 0Pr #°- n 

The equilibrium density and temperature of the gyrokinetic 
species of interest are no and To; the thermal velocity is 
Vth = \/To/m; L is the reference macroscopic scale length 
(i.e., system size), satisfying p/L <C 1 for consistency with 
gyrokinetic ordering. 

In these normalized units, the Maxwellian background 
distribution function is given by Fq = cxp(— (v 2 + 
v 2 )/2)/ (27r) 3 / 2 , and the gyrokinetic equation for the gyroav- 
eraged, perturbed, guiding center density g(R,v ll: v ± ,t) is 
given by 
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where pa — Pth/Pr = Vth/cr = yTo/Tr is the thermal gyro- 
radius pth of the gyrokinetic species normalized to the refer- 
ence gyrorad ius p T . (Our normalization reduces to that used 
in Plunk et al! 35 * 36 ! if we choose T r = To so po = 1, which in 
fact we will do in the 3-D case.) 

The gyrokinetic equation expresses how the guiding centers 
evolve in time due to parallel motion along the magnetic field, 
the gyro-averaged E x B drift across the magnetic field (this 
is the nonlinear term), and parallel electric field acceleration. 
(Note that the slow E x B drift of the guiding center location 
R is different than the rapid gyration velocity v ± of a particle 
around its guiding center.) 

This equation is closed by using the gyrokinetic quasi- 
neutrality equation to determine the electrostatic potential, 
which in Fourier space with these normalized units is given 
by 

0(M) = [ d 3 vJ (k ± p Q v ± )g{k,v ll ,v ± ,t) 
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where 
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f(fc 2 ) = Io(k 2 )e~ k2 is an exponentially-scaled modified 
Bessel function, Io(k 2 ) = Jo(ik 2 ), and r(k) represents the 
shielding by the species that is treated as having a Boltz- 
mann response of some form, the choice of which depends 
on physical situation. If we are treating the ions gyrokinet- 
ically (such as for ion-scale drift waves or Ion Temperature 



Gradient-driven turbulence) and using an adiabatic approxi- 
mation for electrons because of their fast parallel motion rel- 
ative to a typical frequency, k n v te 3> u> (except for modes 
with fe,, = 0), then r(k) = (T r /T e )(l - £ fc|| ), where T e is 
the electron temperature and the discrete Rronecker 5 func- 
tion ensures that the electrons do not respond to zonal modes 
with fc|| = E n = 0. If we are treating electrons gyrokineti- 
cally (such as for small electron scale Electron Temperature 
Gradient-driven turbulence) with an adiabatic approximation 
for ions because k±v t i ^> lj (the k± = mode is not driven 
by any nonlinearities in a periodic domain), then r = T r /Tj. 

Finally, one can also consider a no-response model, t = 
0, which in the 2-D cold-ion limit To — > leads to /? — > 
2-7r/fc 2 , Jo — > 1, and the gyrokinetic equation reduces to 2-D 
hydrodynamics. 

The above difference in zonal flow dynamics for ion vs. 
electron scale fluctuations is responsible for a large enhance- 
ment in zonal flows for ion-scale turbulence, so that zonal 
flows play a key role in the saturation dynamics of ITG 
turbulence^^l and l eads to the Dimits nonlinear shift in the 
critical gradient j^^l j s a i so responsible for a significant re- 
duction in the effect of zonal flows for electron-scale turbu- 
lence, so that they can get to larger amplitude than o ne wo uld 
at first expect from scaling from ion-scale turbulence ! 43 l 44 l 



A. 2D Gyrokinetic absolute equilibria 

For a plasma in a two dimensional (d/dz = 0) cyclic box, 
the collisionless gyrokinetic equation in wavenumber space 
reads 

<9 t <?(k,v)=zx PMpPov)cp(p) -qg(q,v) (5) 

p+q=k 

with the potential (p determined by the quasi-neutrality condi- 
tion 



ip{k) = /3(k) / vdvJ Q (kp v)g(k,v) 



(6) 



where the subscript on v ± has been dropped and the parallel 
velocity v n has been integrated out of the problem. 

The only known rugged (still conserved after mode trun- 
cation) invariants are the "energy" E = (1/2^) J d 2 r[(r + 
(T r /To))ip 2 — (T r /To)ipTip\, and a parameterized set of 
invariants related to the "perturbed-entropy" G(v) — 
(1/2 V) J d 2 Rg 2 (in these equations, V is the volume (area) 
of the integration domain and T is a convolution operator 
in real space given by the Fourier transform of E). (See 
Refs.(4 35 , and 36 ) and references therein for a discussion of 
these conserved quantities and their interpretation.) In Fourier 
space these become E ~ ir ^ k |<,o fc | 2 //3(k) and G(v) = 
J2k u)| 2 /2. As promised in the introductory discus- 
sion, following Leep in what follows we will keep summa- 
tions over only a finite subset K of all possible wavenumbers 
k — the Fourier Galerkin truncation. 45 The Fourier modes in 
the lower half plane are determined by the reality condition, 
g(— k, v) = <?*(k, v), so the state of a system can be uniquely 
specified by the values of the real and imaginary parts of the 



4 



Fourier coefficients g(k, v) for wavenumbers in the upper half 
plane. We will thus consider a further subset K + , defined as 
the modes in K in the upper half plane, which satisfy k y > if 
k x > 0, or k y > if k x < (see also Rrommes and RatrP^j). 
All spectral sums will be expressed in terms of the finite set of 
independent modes in K + , and we denote this summation by 

We can discretize Eq. (|6]l into 



follows: 
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^(k) = /3(k)^ Wi (fc)g(k, Wj ), 



(7) 



where u>j(fc) = miViJo(kp Vi), and is the weight of 
velocity grid point Wj. This discrete form can correspond 
to the case that g(k, v) is uniform on the lattice around 
node i; in general, it is used as a numerical approxima- 
tion for the arbitrary distribution over v as applied in the 
present continuum codes. 46 For a simple midpoint integra- 
tion rule on a grid that extends up to some maximum ve- 
locity v max = wjv, the weight is given by the grid spacing, 
m, = Avi. (More general integration algorithms can also 
be represented in this form. 47 ) With this velocity discretiza- 
tion, there are now N + 1 conserved quantities, given by the 
energy E = Ek 2?r / 3 ( k ) Wig*(k,Vi)wjg(k, vj) and the 

entropy -related quantities (5, = EklffC 4 ' v i)\ 2 - 

So, with the common belief of the applicability of Gibb- 
sian statistical mechanics or Jaynesi31E2l idea of "statistical 
mechanics as a form of statistical inference", we have the dis- 
tribution function ~ exp{— S}, where S is a linear combina- 
tion of conserved quantities, which will be written as 



N 
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Here, the ol; l are the "(inverse) temperature parameters" intro- 
duced as Lagrangian multipliers to form the constant of the 
motion. The Gibbs measure can be shown to be conserved 
by the flow as a generalized Liouville theorem, the incom- 
pressibility of the flow of the phase points in the hyperplane 
spanned by the real and imaginary parts of the Fourier modes. 9 
Note that, importantly, conservation laws and the Liouville 
theorem are inherited, which, with some more assumptions 
(such as ergodicity) makes the discrete system possible to pro- 
duce a Gibbs ensemble. A pedagogical illustration on the 
Gibbs canonical distribution for this system can be found in 
Appendix [A] 

As this is a multivariate Gaussian distribution, one can 
numerically invert the matrix of the quadratic form S = 
Ei,j [$ij a i+ a o'2 7T f3{k)'Wi(k)wj(k)]g(k, Vi)g*(k, vj), but, ac- 
tually, using the Sherman-Morrison formula (see Appendix 
[B| , one can write down the N x N covariances 



c iij (k) = (g*(k,v i )g(k,v j )}/2 
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The spectral energy density D(k) then can be calculated as 
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The isotropic energy spectrum then is E(k) ^ 
spectral density of the perturbed entropy is 
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kD(k). The 
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which expresses the nearly equipartition on Fourier modes 
when the second term in the brackets has negligible 
contribution!^! (Note that D(k) and Gj(k) depend on 
wavenumber not only through /3(k) but also through Wi — 

Wi{k).) 

The temperature parameters are determined by the values 
of the invariants Gi and E. Note that these give the invari- 
ants as nonlinear functions of the temperature parameters ao 
and cti, so this requires numerical solution or further analytic 
approximations to invert. 



Discussion concerning 2+1 D gyrokinetic spectra 

To gain insight into the behavior of these spectra, we will 
explore them in various simplifying limits, such as the small 
and large-fc limits and the cold-ion Hasegawa-Mima limit. We 
will also plot example spectra for particular values of the ao 
and on parameters. 

The gyroaveraging by plasma particles through Eq. ([T]i (fi- 
nite Larmor radius effects) introduces several special func- 
tions whose asymptotic behaviors help shape the energy spec- 
trum. For convenience, let's write here the asymptotic recipes 
for these special functions: 

fc^O, MkpoVi) « 1 - k 2 p 2 vf/4, f (fcVo) ~ 1 - £ 2 Po 



k -> oo, J^(kp Vi) 



2cos 2 (/c«i-7r/2) 
{nkvi) 



T(k 2 pl) 



l 

/2-xkpo 



The simplest limit to consider first is po = \JTq/T t — > 0, 
i.e., the cold-ion limit considered by Hasegawa and Mima in 
their study of drift- wave turbulence 51 (see Appendix [c] for 
more details). The Hasegawa-Mima equation coincides with 
the Charney equation for geophysical flows 52 formulated ear- 
lier, so it is also called the Charney-Hasegawa-Mima (CHM) 
equation. 

In this limit, Jo(kpov) —> 1, and iu,(fc) = 
TniViJo{kpQv) — > rriiVi so the factor of E; u; f(^)/ Q i ln 
Eq. ( fT~0| > becomes independent of k and can just be taken as 
a constant, which we will define as l/(2n 2 o). In the expres- 
sion for /3(k) in Eq. (Wj), we use f (/c 2 p§) — > 1 — k 2 p\ and find 
/3(k) — > 2tt/(t + k 2 ). For comparison with Hasegawa and 
Mima, we set the reference temperature to the electron tem- 
perature, T r — T e , (the electrons are the adiabatic species for 
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ion-scale drift waves) and neglect the Sk\\ factor and associ- 
ated zonal flow effects 53 in the expression for r so that r = 1. 
Note that our normalized k = k ± = k ±>p p s , where k ± p is the 
physical perpendicular wavenumber, and p s — ^T e /nii/ri ci 
is the ion sound radius. (Alternatively, one can consider these 
as equations for electron-scale turbulence where the role of 
ions and electrons is reversed: the ions are adiabatic and a cold 
electron limit is used, in which case the normalizing length is 
an "electron sound radius", p se = y/Ti/m e /tt ce .) 

The result is that the isotropic energy spectrum given by 
Eq. ( p~0] > reduces in the cold ion limit to 

k 



E{k) oc kD(k) oc 



a(l 



2«o' 



(12) 



where a and ao are coefficients that are determined by the 
values of the invariants E and G,;. Note that this is of the 
same form of a 2-parameter family of spectra as in Charney- 
Hasegawa-Mima (CHM) or 2-D Euler absolute equilibrium, 
E(k) oc k/(a C HM + PcHuk 2 ) (Appendix^. 

A remarkable feature of this type of spectrum is that if 
ot-CHM or Pchm are negative, corresponding to a negative 
temperature, then the denominator has the opportunity of 
tending to zero leading to the energy condensation at the low- 
est or highest wave numbers. Energy condensation to the 
lowest modes (one example is shown in Fig. [T| would indi- 
cate an inverse cascade of energy following the argument of 
Kraichnan for the 2D Euler equation 27 (see also Hasegawa and 
Mima 51 , and, Fyfe and Montgomery. 54 ) It then seems that the 
inverse cascade of energy in 2D gyrokinetics could be quite 
a generic feature: Recent r elevant theoretical arguments and 
numerical simulation results 4 -' 3 -^ 5 - 5 -' are consistent with this. 

Instead of the cold-ion limit, we now consider the more 
general case of warm ions (for ion-scale turbulence, or warm 
electrons for electron-scale turbulence), and for simplicity let 
us take T e = Tj = T r (so that po — 1) and neglect zonal flow 
effects (so that r = 1). In the limit k 1, we have /3(k) — > 7T, 
and the magnitude of Wi(k) 2 oc Jy(kvi) will be bounded 
by C/(kVi) for some constant C. Assuming positive on for 
i > 0, we find that the denominator in Eq. ( [T0| approaches 
1 for large k, so D(k) ~ 1/fc and E{k) ~ kD(k) ~ k°. 
This could give a larger tail for gyrokinetics than for CHM, 
which has EcHAi(k) ~ 1/fc (for flcHM > 0). In this 
same k ^> 1 limit, Eq. (Hi simplifies to Gi(k) = l/(2aj), 



which corresponds to equipartition of the generalized entropy. 
We plot the spectra over the reachable wave number regime 
with given temperatures just to sketch the physical picture 
without even bothering to accurately calculate the realizable 
wave number bounds but only with estimations sufficient to 
directly illustrate the problem. (Another way to think about 
the problem and do the corresponding plots is to take the 
bounds of wavenumbers to be prescribed and then realizable 
temperature parameters are determined accordingly. Detailed 
computations and illustrations of more example spectra with 
possible physical discussions are given in Appendix [D] for 
those who are interested.) Actually, much is already known 
from the knowledge of absolute equilibria of 2D EuleP and 
Hasegawa-Mima 51 , though it may still be helpful to give some 
general physical picture, especially the finite Larmor radius 
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FIG. 1. (Color online) Example spectra for various values of ao, 
with a 4 = 10 3 exp{u l 2 /2} for i > 0, and p = 1CT 1 , r = 1. 
Negative ao state can occur that correspond to condensation of most 
of the energy into the longest wavelength modes. 



effects, with some example spectra as shown in Fig{T] The on 
for i > were set by a^ 1 = ICC 3 cxp(— vf/2). We take 
N = 40 with Vi homogeneously collocated between and 
V = 3, and, K = {k|l ^ k x ,y ^ 150} for positive tem- 
peratures cases while K = {k|4 ^ k XjV ^ 150} for a neg- 
ative temperature case. Some details, including the values of 
rrii (= 1 here), are not important, and reasonable changing 
of them (as a re-normalization of the variables) won't affect 
our results. Here, as E(k) ~ fcD(k), the low-A; equiparti- 
tion range have E(k) ~ k and E(k) ~ k° for large fc, both 
of which can be easily checked with the asymptotic recipe of 
the special functions given in the beginning of this subsection. 
The transition in between represents the FLR effects, the de- 
tails of which also depend on the details of the temperature 
parameters. Other temperature parameter values may change 
the k ranges where such asymptotic behaviors can be realized; 
or, if the truncation wave numbers in the computations k m i n 
and k max (which may be relevant to some characteristic, such 
as the collisional, scales in real physical systems) were not 
chosen properly, we would not be able to reach such asymp- 
totic behaviors. 

A negative value of ao corresponds not only to an enhance- 
ment of energy at larger spatial scales (low k) but also to an 
enhancement of fluctuations at larger velocity scales, as given 
by Eq. (|9j. If ao — 0, then this equation indicates that differ- 
ent velocity grid points i ^ j are uncorrected, while making 
a more negative will increase the correlation length in the 
velocity, particularly at low fc. If collisions are included, they 
will cause dissipation at both small velocity scales and small 
spatial scales (through FLR effects corresponding to classical 
diffusion),^ so the inverse cascade found here in 2D will tend 
to reduce both forms of dissipation. 

We end up this sub-section by remarking that, since the 
CHM limit absolute equilibrium statistics has already been 
verified by numerical experiment] 5 -^ our theoretical calcula- 
tion for gyrokinetic is then also, to some degree, endorsed. 
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B. 3+2D Gyrokinetic absolute equilibria 



for S — jS — j(E + W g o) gives 



The mathematical treatment for the calculation of Galerkin 
truncation absolute equilibrium is basically the same for 3+2D 
and 2+ ID cases. However, some brief remarks about the sys- 
tem and the conserved quantity, followed by some technical 
details in the calculation are still necessary. 

In the full gyrokinetic equation with three spatial and two 
velocity dimensions, notice that two linear terms, from paral- 
lel motion along the magnetic field and from the parallel elec- 
tric field acceleration, are simply added to the equation for the 
(2 + 1)-D case, without changing anything about the nonlin- 
ear term. While the gyroaveraged E cross B drift conserves E 
and G separately, the parallel motion makes them talk to each 
other and combines them into another conserved quantity, the 
generalized energy W = E + W g o (see Refs.dH[33][36] and 
l57l l and references therein for further discussions of this quan- 
tity): 



W 



err 

w 



[(1 + T)lfi 2 - tpTif] + / d 3 V 



g 



2V F 



(13) 



(For simplicity in all of our 3+2D work, we will set the refer- 
ence temperature T r used in normalizations to the temperature 
of the kinetic species To, so po = 1.) The Fourier Galerkin 
truncated form of the generalized energy is 



S 



|^?(k)| 2 + 7r // v±dv±dv 



l3(k,v) 



We will discretize velocity space in a way that makes it easy 
to reduce the previous 2D spatial + ID velocity results to the 
3D spatial + 2D velocity results here. Specifically, we will 
discretize velocity integrals as 



d vg(v) = 2tt / dv ± / dv ll v ± g(v ll ,v ± ) 

JO J-oc 
N 

m 2K^miV ±ti g(vi), (14) 



where i now indexes over all points v, = (v ± i,v^) in the 
2D velocity space grid. For a logically rectangular mesh in 
(v±, there would be a total of N — N V± N V]] grid points, 
with N v± points in the perpendicular velocity direction and 
N v „ in the parallel velocity direction. (As in the ID ve- 
locity case, for a simple midpoint integration algorithm, rrii 
is the weight of the i'th velocity cell, rrii — Av^At!^*, 
while more generally the weights m t and grid point locations 
(v±.i, I'll,,) can be chosen to give high-order Gaussian quadra- 
ture.) The 2D velocity generalization of Eq. |7]), the dis- 
cretized quasineutrality equation to determine the potential, 
now reads <£>(k) = /3(k) J^i VJ i(k±)g(k, Vj), where Wi(k ± ) — 
miVx,iJo(k x v ±ti ). 

We then can calculate the absolute equilibria following the 
same procedure as in the 2D case, but now only one inverse 
temperature parameter 7 shows up in the canonical distribu- 
tion ~ exp{— 75}. Using the above velocity discretization 



S = 72^ k 27r/3(k)E ij w; i (?*(k,v i )w i 5(k,v J -) 

+^*L 1 27rm I v ± ^ k \g(K^)\ 2 /F (v t ) 

We note from this that a negative temperature is not re- 
alizable any more. Compar ing this expression for the 3D 
S with the 2D result in Eq. ( Al 1, we see that they become 
identical if we make the substitutions uq = 7 and Q4 = 
2n'ymiV ±i i/Fo(\i). All of the 2-D results thus generalize to 
the 3-D case with these variable substitutions. For example, 
the electrostatic component of the spectral energy density in 
Eq. ( p~0] > becomes 



D(k) 



1 

2^ 



/3(k) J2i m i v ±,i F o('"i) J o( k ± v ±,i) 

1 + /?( k ) E 4 miVxjFo^JfiftxVxj) 



the small lattice size limipil, where 

use 2TTj2 l m i v ±i iF (v i )J^(k ± Vx,i) 



(15) 



we 



In 

can 

J d 3 vF (vi)Jo(k ± Vx,i) — r (fc 2 ), the electrostatic po- 
tential spectral density becomes 



r (fci) 



7 (r + l-r (fci))(T + l) 



(16) 



The shape of this spectrum is consistent with the discrete- 
particle thermal noise spectrum for gyrokinetic PIC codes cal- 
culated by one of us previously, as given in Eq. (5) of Ref. d34i l. 
which reduces to the above result in the limit where numeri- 
cal details such as spatial filtering and finite differencing^ are 
ignored by setting 5c(k) = 1 and G?n(k) = 1, and by taking 
the r = 1 limit in our expression. (The thermal spectrum in 
Re f. (1341 was calculated for the case of one gyrokinetic species 
and one adiabatic species, as also assumed in the present pa- 
per, and also accounted for various numerical factors as used 
in typical PIC codes. The first calculations of the discrete- 
particle thermal noise spectrum for gyrokinetic particle codes 
are in Refs.dT6landl59il and were for the case where all species 
were treated gyrokinetically.) 

Ref. d34b found good agreement between this analytic ther- 
mal spectrum and the fluctuation spectrum in a PIC code in 
a noise-dominated regime, providing support for the calcula- 
tion done here. Readers interested in a discussion of noise in 
numerical schemes are referred to Appendix El 



Relevance to 3D plasma turbulence 

There are several interesting features of the 3D spectrum in 
Eq. ( 16 1. Note that it is independent of fc||, i.e., the equilib- 

so presum- 



rium spectrum corresponds to equipartition in k n 
ably the nonlinear dynamics of a turbulent system should tend 
to drive cascades to high fcy. (Gyrokinetics assumes fc y <C k ± , 
so there is a limit to how far this spectrum can extend within 
this model.) Also note that even with the finite-Larmor ra- 
dius averaging in gyrokinetics, the electrostatic potential spec- 
trum falls relatively slowly at high k ± since Tq ~ C/k ± , so 



7 



the electrostatic energy spectrum is flat at high wave number, 
E<p(k±) oc£^(k)| 2 

For ion-scale non-zonal flows with adiabatic electrons, the 
long-wavelength limit of Eq. (16i is (l^iJ 2 ) = 1/(27) ( se t~ 
ting t = T r /T e = Ti/T e = 1 for simplicity). But for 
zonal flows, which have k y = k z = fc B = and thus have 
t = (see the discussion after Eq. ([4])), the long wavelength 
limit is (\<Pzf k| 2 ) = V il^x) (h ere the "ZF" subscript refers 
to the zonal flow component of the potential), so the ampli- 
tude of long-wavelength zonal flows is enhanced relative to 
other nearby modes by a factor of ~ l/k x . (While the re- 
sulting zonal potential blows up as k x —> 0, the shearing rate 
oc dvy/dx oc d 2 (p ZF /dx 2 oc k x ip ZFk oc k x remains well- 
behaved.) However, one of us, GWH, tends to believe that 
this enhancement in the 3-D statistical equilibrium is interest- 
ing but by itself is probably not enough to explain the observed 
importance of zonal flows in ITG turbulence, since there are 
very few zonal modes compared to the many other modes with 
k y ^ or fcy 7^ 0. The importance of ITG zonal flows is prob- 
ably due to other effects, such as the way in which the lack of 
adiabatic electron response causes an enhancement of the sec- 
ondary instabilities^ (or related parametric instabilities) that 
drive zonal flows. 

However, much stronger enhancement of zonal flows can 



exist in the 2-D absolute equilibrium of Eq. ( 10 1 where nega- 
tive ceo can strongly enhance modes with r = 0. The mech- 
anism for this enhancement is related in a way to the en- 
hancement of zonal flows in secondary/parametric instabilties. 
This 2-D equilibrium effect might be related to the enhance- 
ment of zonal flows in an actual turbulent plasma, if there 
are regions of the turbulent spectrum where the parallel dy- 
namics is slow compared to the nonlinear decorrelation rate 
fcyVt -C Aojatl ~ k ± VExB and so act in a quasi-2D man- 
ner. However, there will also be competition from 3D effects, 
which limits the inverse cascade and tends to push the spec- 
trum towards equipartition in fey. 

The 2D and 3D gyrokinetic absolute equilibrium results 
may also provide insight into other aspects of driven non- 
equilibrium gyrokinetic turbulence, such as the directions of 
turbulent cascades in (k ]u k ± ). The inverse cascade found in 
2D may imply that in regions of a turbulent spectrum where 
fe||Uf <C k ± VExB, then the interactions may be quasi-2D and 
undergo an inverse cascade to smaller k±, simultaneously with 
a cascade to higher fcy (towards equipartition in k n ), until the 
parallel dynamics becomes competitive with nonlinear terms, 
k\\V t ~ k ± VExB- At this point it might then switch to a 
forward cascade to higher wavenumber, but along a path in 
(fey, k ± ) space such that k n vt ~ k ± VExB- Thus this supports 
the critical balance hypothesis suggested for gyrokinetic tur- 
bulence in Refs.(4 and 33), that the turbulence will primarily 
cascade along a path in wave number space that has parallel 
linear time scales comparable to perpendicular nonlinear time 
scales, similar to critical balance ideas in astrophysical Alfven 
turbulence in Refs.(l5T1andl62l. Further analysis of gyrokinetic 
statistical equilibria may lead to more specific insights. 

There are other more subtle physics, such as the bottleneck 
and its associated weakening of intermittency growth issues,^! 
as proposed to be explained as partial thermalization by Frisch 



et alP For example, as the Fourier transform is linear, the 
physical-space field of the Fourier Galerkin truncated absolute 
equilibria would also be Gaussian, whose residual may result 
in a resistance in the departure from Gaussian (intermittency) 
for the turbulence fluctuationsP^Before examining the details 
of collision and wave-particle interaction mechanisms, so far 
we unfortunately are not able to say anything more on this for 
the plasma turbulence. Nevertheless, such considerations em- 
phasize the importance of implementing the correct collision 
operators (which is necessary in many physical situations) and 
in interpreting the numerical data. 



III. CONCLUSION AND FURTHER REMARKS 

Here we have extended previous work on statistical equi- 
libria of 2D and 3D hydrodynamics and MHD to the case of 
higher-dimensional gyrokinetics. Previous work in hydrody- 
namics found that there was a profound difference between 
2D and 3D, because the existence of 2 invariants in 2D lead to 
the existence of negative temperature equilibrium states with 
most of the energy condensing into the longest wavelengths in 
the system (related to the inverse energy cascade in 2D turbu- 
lence), while in 3D there was only a single invariant resulting 
in energy equipartition among Fourier modes (related to the 
forward cascade of energy to small scales in 3D turbulence). 

For gyrokinetics in the limit of 2 spatial and 1 velocity di- 
mension (2+ ID), we have worked out the Gibbs equilibrium 
in the presence of N + 1 invariants (where N is the number 
of velocity grid points) and find that, like 2D hydrodynam- 
ics, this can also exhibit negative temperature states where 
much of the energy condenses to the longest wavelengths 
in the system. For a range of typical parameters explored 
so far, 2+ ID gyrokinetics exhibits a very strong inverse cas- 
cade. At high k ± , the 2D gyrokinetic energy spectrum has an 
asymptotically-flat tail, E(k ± ) ~ which is enhanced rel- 
ative to the high k ± limit of Hasegawa-Mima's thermal spec- 
trum, E(k ± ) ~ l/k x . The amplitude of this tail in gyrokinet- 
ics is found to depend sensitively on the ratio of Gi to energy. 

We also calculated the statistical absolute spectrum for 
Fourier-truncated gyrokinetics in the full 3 spatial and 2 ve- 
locity dimensions, and found that the result was equivalent 
to earlier thermal noise spectra calculated for particle-in-cell 
gyrokinetics, indicating that the random phase and amplitude 
of shielded Fourier components of the distribution function in 
a continuum representation is related to the random position 
and weights of shielded particles in the Klimontovich repre- 
sentation of a PIC code. The resulting 3-D gyrokinetic spec- 
trum corresponds to equipartition in k n , and even with all of 
the finite-Larmor radius averaging in gyrokinetics, the elec- 
trostatic potential spectrum only falls relatively slowly at high 
k ± , so E(p(k±) oc k ± \ip(k)\ 2 ~ k°. 

As described in the introduction, statistical equilibria spec- 
tra as calculated here have several useful purposes. In partic- 
ular, they provide an analytic nonlinear test for benchmarking 
of gyrokinetic codes, which could be pursued in future work. 
They may also provide insights into certain aspects of the non- 
linear dynamics in driven, non-equilibrium gyrokinetic turbu- 
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lence simulations. For example, in regions of the turbulent 
spectrum where the parallel linear dynamics is slow compared 
to the nonlinear decorrelation rate, k n v t -C k ± VE X B, then the 
interactions may behave in a quasi-2D behavior, which can 
cause an inverse cascade to smaller k ± in general, and in par- 
ticular can strongly enhance the ITG zonal flows because of 
the lack of adiabatic electron shielding for ITG zonal flows. 
But these may be offset by the tendency towards equipartition 
of the spectrum in k n , so that eventually k ]f v t ~ k ± VExB and 
parallel linear dynamics becomes competitive with nonlinear 
perpendicular dynamics. In this region of wavenumbers, the 
turbulent cascade would then switch to a forward cascade to 
higher |k|, along a path where parallel and perpendicular dy- 
namics remain comparable and so stay full 3D, consistent with 
the critical balance hypothesis for gyrokinetic turbulence sug- 
gested in Refs.(4] and [33l l. There are various directions in 
which the present work could be extended in the future that 
may further help in understanding plasma behavior in actual 
experiments, such as extensions to include a kinetic treatment 
of all particle species, electromagnetic fluctuations, and the 
effects of magnetic curvature and grad-i? drifts in toroidal ge- 
ometry. 
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Appendix A: Pedagogical illustration of the calculation of 
the canonical ensemble 

The line of reasoning presented by T.-D. Lea2l regarding 
how to apply a statistical mechanics approach to hydrodynam- 
ics and MHD can be straightforwardly extended to the higher 
dimensional gyrokinetic case considered here. We briefly 
summarize that line of reasoning here, which also serves to 
explain the notation that we use. Consider a system governed 
by Eqs. |5]l and |7} (with Eq. |5} evaluated at the same velocity 
grid points as used in Eq.(|7])). The state of a system at a partic- 
ular time can be specified by a vector g in an extended phase 
space of dimension N^N, where is the number of Fourier 
modes and N is the number of velocity grid points. The el- 
ements of g are g(k, Vi), the complex amplitude of Fourier 
modes k € K + at velocities vu where K + is the set of inde- 
pendent modes k in the truncation K that are in the upper half 
plane. One can consider an ensemble of many such systems, 
and define the function V(g, t) that gives the probability of a 
system being in state g at time t. For continuous dynamics, 



this satisfies a conservation law d{P + d s ■ (gP) = where 
an over-dot is used to denote a time derivative so g is given 
by Eq. Q. A generalized Liouville theorem holds for these 
equations, i.e., the flow in this extended phase space is incom- 
pressible, <9 g • g = 0, because for a given value of k, the right 
hand side of Eq. |5]l vanishes if p = ±k (because q = k p 
means pxq vanishes), and thus also vanishes if q = ±k. (In 
other words, the rate of change <jt(k) at any instant in time de- 
pends only on the amplitude of other modes g(p) with p ^ k.) 

Since a Liouville theorem holds, standard results and as- 
sumptions from statistical mechanics can be applied to these 
equations. A generalized Liouville equation holds, d{P + 
g • d g P = 0, i.e., the probability V(g(t),t) is constant on a 
moving trajectory in this extended phase-space. Looking for 
a time-independent statistical steady state, we take an equal 
probability for all points along a trajectory's path. Assuming 
that the dynamics are sufficiently mixing and an ergodic hy- 
pothesis holds, so that a trajectory samples all possible points 
on a hyper-surface in phase-space constrained only by the 
invariants, leads to the micro-canonical ensemble given by 
V = CS{E - B )niI 1 <5(G i - Guy), where E = E(g) and 
Gi(g) are the previously given expressions for the energy and 
entropy invariants, which are functions of g, Eq and G,o are 
the values of those invariants (set by initial conditions), and 
Tlf—x indicates repeated multiplication over all possible ve- 
locity points i. 

As is well known, for systems with a large number of de- 
grees of freedom, many features of a micro-canonical en- 
semble are often well-approximated by a Gibbs canonical 
ensemble, V = Z _1 exp(— S) where S is a linear combi- 
nation of conserved quantities, which in this case is S = 
a^E + ^ i otiGi, ao and the N values of oti are the "(in- 
verse) temperature parameters", and Z is a normalization co- 
efficient such that J dgP(g) = 1. One wa> l 48 l 49 l to derive 
this is to choose V to maximize the Liouville phase-space en- 
tropy Sl = — J dg'P(g) log("P(g)) (i.e., choose V to be as 
uniformly distributed as possible) subject only to constraints 
on the average values of the invariants (this leads to the La- 
grange multipliers o.q and a,; in the canonical ensemble). For 
example, the constraint on the ensemble-averaged value of the 
energy is E = (E) = J dgV(g)E(g). 

[Note that the fact that a Liouville theorem is satisfied is an 
important part of justifying the maximum-entropy approach 
of the previous paragraph, as it means that the probability dis- 
tribution V can be constant along trajectories in these coordi- 
nates, which is not necessarily true in other coordinates. For 
example, something that is uniformly distributed in x is not 
uniform in a; 3 .] 

Inserting the expressions for the energy and entropy invari- 
ants into the the expression for S gives 

S = ao^ k 2^(k)^r(k,^)^(fc)^(fc)g(k,^) 

N 

+ E^E k l5(k,^)| 2 (Al) 

i=l 

With a little rearrangement, the Gibbs canonical distribu- 
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tion becomes 



we have 



7>«exp(--]T k g*(k)-M(k)-g(k; 



(A2) 



This is of the form of a multivariate Gaussian distribution, 
where the elements of the k-dependent, N x N matrix M are 
given by — 5ij2cti + a 'iTr(3(k)wi(k)wj(k), and g(k) 
is the N -dimensional vector of the velocity-indexed values of 
the complex amplitudes g(k, vi). 

Expressing g in terms of its real and imaginary parts, 
g(k) = gj?. (k) + igi (k), note that the sum over wavenumbers 
in Eq. (A2i can be written as 5Z k g.n:(k) • M(k) • gij(k) + 
^ k g/(k) • M(k) • g/(k) since M is real, so the real and 
imaginary parts of g are uncorrelated and have the same 
co-variance, (g R (k,Vi)g R {k,Vj)) = (gi(k,Vi)gx(k,Vj)) = 
Cjj(k), where cy are the elements of the co-variance matrix 
C(k) given by the inverse of M, i.e., C(k) = M _1 (k). 



Appendix B: Calculating the covariance matrix using 
Sherman-Morrison formula 



In principle, once that M(k) in Eq. (A2i is known, one 
can calculate the co-variance matrix C(k) = M _1 (k), and 
one can then calculate various statistical properties of interest, 
such as the energy spectrum of fluctuations. However, M is in 
general a dense matrix, so at first it looks like this may require 
a numerical treatment to invert it. To make analytic progress, 
one can initially consider the limit ao = 0, in which case 
M is diagonal and easily invertible. One can then do a matrix 
series expansion for small a and discover that it is possible to 
sum the result to all orders in ao because of the special tensor 
product form of the coefficient of ao in M. This turns out to 
be a special case of the general Sherman-Morrison formula. 

In linear algebra, suppose A is an invertible square matrix 
and u, v are vectors and that 1 + v T A^ 1 u ^ 0, then the 
Sherman-Morrison formula reads 64 



A~ 1 uv T A- 1 
1 +v T A- 1 u' 



To derive the covariance matrix C = M _1 in Eq. JfS, 
we use the definition of M given after Eq. ( |A2| i. Thus 
in the Sherman-Morrison formula, the elements of A 

are a^j — 2on§ij, and we can set u\ — ao4ir(3wi 
and Vi — itfj. Then the elements of A~ l uv T A -1 
are Xij — aQ2Trf3{k)wia~ 1 WjOi~ 1 , and v T A~ 1 u — 

a 27r/?(k) J2 t Y,j w i w j a~ 1 5 lj = a 27r/3(k) £\ wf^ 1 - So > 
we have Eq. d5). 



Appendix C: From gyrokinetics to fluids: recovering the 
Charney-Hasegawa-Mima equations 

We first briefly reproduce the derivation of the Charney 
Hasegawa-Mima (CHM) equations from gyrokinetics 
Plunk et alP^ with slight variation: From Eqs. ^ and |6 



d t <p(k) = /3(k) i x P 1^(P) 

p+q=k 

x vdvJ (kp v)J (pp v)g(q,v). (CI) 



In the cold ion limit po — > (as described in Sec. |II A| >, 
the first kind zeroth order Bessel functions reduce to unity, 
and /3(k) to 2tt/(t + k 2 ), and then, with substitution of the 
quasi-neutrality Eq. |6]l in the second line of (CI I, gyrokinet- 
ics reduces to CHM. In physical space, it reads 



d t {r - V 2 )ip = zxV^ V(VV) 



(C2) 



This is the inviscid (the collision operator in Plunk et alP^also 
vanishes after integration over velocity by particle conserva- 
tion) CHM equation. The scale to which gradients were nor- 
malized in these equations corresponds to the Rossby defor- 
mation radius in quasi-geostrophic turbulence, or to the sound 
Larmor radius, p s , in a plasma. We have left a r dependence 
in these equations for generality, as the two-dimensional Euler 
equation can be obtained in the case r = (the no-response 
model). 

There are two invariants of the CHM/Euler equation that are 
relevant to our discussion, referred to as energy and enstrophy 
(although their physical interpretation depends on the specific 
scale of interest): 



E, 



CHM 



J CHM 



cPr 
V 



[V + |V^| 2 ] (C3a) 
[t\V<p\ 2 + (V 2 <p) 2 }. (C3b) 



This leads to absolute equilibrium energy spectra of the form 
E(k) oc kD(k) oc k/(a C HM + PcHNik 2 ), where occhm and 
Pchm will be determined by the values of these two invari- 
ants, via Echm = 2E kJ D(fc) and Z C hm = ^E k k 2 D{k). 
The first invariant, Echm, is formally the reduced energy, E, 
of gyrokinetics in the cold-ion limit. The enstrophy, Zqhm, 
however is new and deserves further inspection of its origin. 

The CHM equation can be written asdn/dt — zxV<p- Vn, 
where the potential is determined from the guiding-center 
density n = 2ir J dvvg by inverting (t — V 2 )ip = n. The 
nonlinear term on the RHS of CHM has the property that 
it can be multiplied by either the density n or the poten- 
tial <p and then will vanish when integrated over all space. 
This leads to the two standard energy and enstrophy invari- 
ants used for the CHM equations. However, in gyrokinet- 
ics where we keep a finite, non-zero temperature, the veloc- 
ity and wavenumber dependence in the Bessel functions in 



the second line of Eq. (CI I introduces extra non-local posi- 
tion and velocity scale interactions that mean that Zchm is 
no longer conserved. Due to these FLR effects, gyrokinet- 
ics instead has a set of invariants that hold at each velocity, 
G(v) oc J d 2 Rg 2 (R, v), while CHM had an additional invari- 
ant proportional to J d 2 R(J dvvg) 2 . This new CHM invari- 
ant is not representable as a combination of the gyrokinetic 
G(v) and E invariants. One way to think of this is to note 
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that CHM depends only on the velocity integral of g through 
n = 2ir J dvvg, so the CHM dynamics are independent of 
any details of the velocity structure of g, thus allowing an ad- 
ditional invariant that is not present in gyrokinetics because of 
its FLR effects. 



Appendix D: Example 2D Spectra for Specified Initial 
Condition^ 



Here we consider a numerical gedankenexperiment, in 
which a gyrokinetic code is operated in the 2+ ID limit and 
is initialized with perturbations concentrated near some initial 
wavenumber fc but with no other forcing. Those perturba- 
tions will then interact nonlinearly and scatter energy to other 
wave numbers, while preserving certain invariants of the mo- 
tion. Presumably the spectrum will eventually reach a statis- 
tical steady state, and here we make plots of the energy spec- 
tra expected from canonical equilibria corresponding to some 
sample initial conditions. This helps provide further insight 
into the nature of these equilibria. 

Before making these plots, we first consider some of the 
properties of spectra and the relationship between the invari- 
ants and the parameters in more detail. For positive ao 
and a>i (i > 0), then the factor in brackets in Eq. (Hi is 
close to unity for all wavenumbers, and one can sum over 
all wavenumbers to find Gi = N^/(2ai) (where Nk « 
^(kmax/kmin) 2 is the number of Fourier modes), which can 
be used to determine on in terms of the conserved G\. Eq.( 10 1 
can be summed over all wavenumbers to determine the total 
energy and then determine a . For fixed positive values of 
on, the energy is a monotonically decreasing function of ao, 
so if the energy is sufficiently large (for given values of the 
Gi), then ao must go negative to produce a "negative temper- 
ature" state. If we perturbatively use a.; oc 1/Gi to evaluate 
the energy spectrum, and assume a Maxwellian velocity dis- 
tribution for the fluctuations so Gi = J d 2 Rg 2 (TL, v)/2 oc 



i w i ' 



exp(— v 2 ), then the commonly occurring factor 

v ?Jo(kvi) exp(— v 2 ) is a monotonically decreasing func- 
tion of k, as is /3(k), so if ao goes negative in the denominator 
of Eq. ( fT0] i, it will preferentially enhance the energy in the 
low-fc part of the spectrum. (If the denominator of Eq. ( 10 1 
gets sufficiently close to zero for some wavenumbers, then the 
factor in square brackets in Eq. (JTTJ could differ from unity 
and alter the relationship between Gi and assumed here.) 
Note that the realizability constraint that the energy spectrum 
be non-negative in this case means that the limiting value of 
ao for this set of a/s is au m — — [27r/3(k) J2i wf(k)aY]~ 1 
evaluated at k = k min . (If the enhancement of ion-scale zonal 
flows due to the lack of electron response is accounted for, 
then this would strongly increase the value of /3 for the zonal 
modes, reduce the magnitude of the limiting value of ao, and 
strongly enhance the amplitude of zonal flows.) 

Returning to the numerical gedankenexperiment, consider 
an initial perturbation of the form 

g{R,v,t = 0) = cos(k R y )— — J (k v) (Dl) 

Z7T 



(here we set po = 1 for simplicity), as a model that has some 
characteristics of the drive by drift-wave types of instabilities. 
This initial condition models what happens if a linear source 
term —vexb ■ Vi<o (representing instabilities that drive drift- 
wave gyrokinetic turbulence) had been turned on in the gy- 
rokinetic equation for a time of order L/c r in the presence of 
a background density gradient VF = —xF /L, where the 
eddy has a bi-normal wavenumber k y — ko- (In an actual 
code, a small amount of energy must initially be put in other 
Fourier modes as well, because a single Fourier mode does 
not interact with itself nonlinearly.) The energy and entropy 
invariants corresponding to this initial condition are E = 
/3(fc )f 2 (fc 2 )/(8^) and G % = exp(-w 2 ) J 2 (fc ^)/(167r 2 ). 

Given the specified values of the energy and entropy in- 
variants, it is not analytically easy in general to invert the 
equations to determine the corresponding temperature pa- 
rameters ao and a^, because the energy and entropy are 
nonlinear functions of the temperature parameters, as dis- 
cussed after Eq. (111. That is, the entropy invariants are 
related to the covariance matrix by G; = 2^ k c^ j (k), 
and the energy is related by E = 27r^ k |c/? fc | 2 //3(k) = 



N 



.j(k), where the wavenumber sums 



E k 47r/3(k)E 

are over the independent set K + and Cjj(k) is a nonlinear 
function of the temperature parameters as given by Eq. |9]). 

A code was written to numerically carry out the inversion 
using a nonlinear root solver based on Powell's method and 
Broyden's quasi-Newton algorithm in the minpack software 
package. 66 To aid in finding a root, a variable transformation 
was used for ao to ensure that during the search ao never 
exceeded the lower limit set by realizability constraints that 
the energy spectra be non-negative. The results in this sec- 
tion used a uniformly spaced 2-D wavenumber grid, 67 where 
the set of retained Fourier modes is K = {k| k m i n < |k| < 
k max + Afc/2}, with k min = Ak = 0.1, k max = 10, and 
the number of modes is Nk = 31, 576. A uniformly-spaced 
velocity grid was used with N = 40 points, equally spaced 
from Aw/2 to v max — 3 + Av/ 2, with m, = Av — 3/N. The 
background plasma temperatures were set to T = T e = T r so 
po = 1 and an adiabiatic species response factor of t = 1 for 
simplicity, neglecting possible enhancements of zonal flows. 

Fig. |2]) shows the gyrokinetic equilibrium spectrum that 
results from these initial conditions with ko = 0.3. (The 
wavenumbers in the figures refers to the physical k p , where 
the normalized k = k p p r , and the reference gyroradius is 
set to p s = \jT e jniil^l C i for comparison with the cold- 
ion Hasegawa-Mima drift-wave equations.) This figure also 
shows the spectrum given by the Charney-Hasegawa-Mima 
equations for these same initial conditions. 68 Both the gyroki- 
netic and CHM spectra show strong transfer of energy to large 
scales relative to the initial location of the energy at ko = 0.3, 
though there is more of a tail in CHM case. Both the gyroki- 
netic and CHM equations result in a negative temperature state 
(ao < for the gyrokinetic case 6 -?) with most of the energy 
condensed into the longest wavelengths in the domain. 

Fig. Q is similar to Fig. Q except that the energy is ini- 
tially at a higher wavenumber of fc = 5.0. There are now sig- 
nificant differences between the gyrokinetic and Hasegawa- 
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FIG. 2. (Color online) Spectra for 2+1D gyrokinetics (GK) and 
for the Charney-Hasegawa-Mima (CHM) equations, corresponding 
to the model initial conditions with energy initially at kop a = 0.3. 
Both spectra show a significant transfer of energy to larger scales, 
resulting in a negative temperature state with most of the energy con- 
densed into the longest wavelength in the domain. 



Mima spectra, with gyrokinetics still showing a very strong 
transfer to large scales while the energy remains primarily at 
higher k in the Hasegawa-Mima case. This is because the 
cold-ion Hasegawa-Mima equations have an additional invari- 
ant (see Appendix O, the enstrophy (the mean squared vor- 
ticity), which is not conserved by the general warm-ion gy- 
rokinetic equations because of FLR effects in the Bessel func- 
tions. (This is related to the fact that although the 2+1D gy- 
rokinetic spectrum in Eq. ( 12 1 for the T /T r <C 1 regime has 



FIG. 3. (Color online) Spectra for 2+1D gyrokinetics and Charney- 
Hasegawa-Mima, like Fig. l|2j except the energy is initially at a 
higher wavenumber of kop s = 5.0. 



the same 2-parameter form as the Charney-Hasegawa-Mima 
(CHM) spectrum, the relationship between those 2 parame- 
ters and the invariants is different for gyrokinetics than for 
CHM, because CHM has an additional invariant at To = 
that doesn't exist in gyrokinetics with non-zero To.) The ini- 
tial wavenumber of ko = 5 in Fig. ^ is sufficiently close to 
the truncation wavenumber k max = 10 that there is not much 
room for the enstrophy density oc k 2 E(k) to transfer to higher 
wavenumber, thus inhibiting how much transfer of energy to 
larger scales can occur in the Hasegawa-Mima equations. 

It is possible to increase the size of the tail in the 2+1D gy- 
rokinetic spectrum by increasing the amplitude of the Gi rel- 
ative to the energy, as shown in Fig. Q. Considering a long- 
wavelength initial condition ignoring FLR effects, this can oc- 
cur if a component is added to g(R, v,t = 0) that oscillates 
in velocity so that it makes no contribution to the potential oc 
J dvvg, but does enhance Gi — J d 2 R \g(R, v.i)\ 2 . This can 
model the effects of temperature gradients in the background 
Fq that drives the initial perturbation, or the build up of large 
values of G; in a long turbulence simulation without adequate 
dissipation because of the entropy balance relationships!^^ 
thus leading to bottleneck problems. 29 Note that the depen- 
dence of the tail on the enhancement of Gi is a strongly non- 
linear functionp2lln the limit of very large Gi / E, the spectrum 
will approach equipartition among Fourier modes, E oc k. 
One can also consider how the spectrum depends on the as- 




FIG. 4. (Color online) Spectra for 2+1D gyrokinetics for initial con- 
ditions with the energy at kop s — 0.3 like Fig. Q, but with the 
value of the entropy invariants Gi enhanced relative to the model 
initial conditions by factors of 1 x , 4 x , lOx, and 500 x . 



sumed velocity grid spacing Av. From numerical results, con- 
firmed by analytic scalings, one finds that as Av goes to zero, 
with fixed values of E and Gi, that ao oc — 1/ Av goes to neg- 
ative infinity (while on — > constant for i > 0), corresponding 
to a negative temperature state with all of the energy in the 
lowest k mode, so E{k) — for all k > k m i n . 



Appendix E: Thermal noise spectra in numerical schemes 



In the 3+2D results in Sec. dllBb, we worked out the elec 



trostatic spectrum, Eq. ( 16 1 and showed that the shape agrees 
with earlier results for the spectrum in a PIC code. Here we 
show that the magnitude agrees as well, with the proper re- 
lation between certain quantities in a continuum code and a 
PIC code. Begin by defining a weighted mean-square average 
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of the distribution function g 2 = J d 3 R J d 3 v(g 2 ) /(F V) (an 
overbar is used here to indicate a combined velocity space av- 
erage and an ensemble/volume average, to be distinguished 
from angle brackets that indicate an ensemble average). This 
uses the same velocity weighting as found in the W g o compo- 
nent of the generalized energy in Eq. ( fj"3| ). After discretiza- 
tion, this becomes g 2 = 2^ k ^ ™zCi,i(k)/Fo(v;). Us- 
ing Eq. ( [TT] ) evaluated with the coefficients given just before 
Eq. (15 1 for the 3+2D case, and using the same approxima- 



tions as used just before Eq. ( 16 1 (where velocity summations 
are approximated by integrals assuming a well-resolved ve- 
locity limit), one can show that 



constant set by initial conditions, but in the case of turbu- 
lence driven by a background gra dient, g 2 will increase in time 
due to an entropy balance relatiorP^ISZIunless dissipation is ac- 
counted for.) Most continuum codes avoid such problems by 
employing either physical collisions or numerical dissipation 
such as high order upwinding or hyperdiffusion, though work 
on improved subgrid models might be able to help optimize 
the performance by reducing resolution requirements. While 
6f PIC simulations can formally work correctly for a given 
simulation time period if enough particles are used, eventu- 
ally the noise can grow in time to become a problem. PIC 
codes can avoid this issue and/or reduce the particle resolu- 
tion requirements by employing weight-resetting methods,^ 
which essentially provide some numerical diffusion to limit 
the growth of the weights. 



for N 3> 1 (recall that ^ k i s defined as a sum over the modes 
in the upper half plane, so the total number of Fourier modes is 
2^k = Nft). We thus find that the thermal noise spectrum in 
Eq. (5) of Ref.|34]for §f PIC codes (using a weighted-particle 
Klimontovich representation for the distribution function) is 
identical to the thermal spectrum calculated here for a contin- 
uum code using a spectral representation for the distribution 
function, with the identification of 1/7 = g 2 / (NNk ) in a con- 
tinuum code with w 2 /N p in a PIC code. So the total number 
of particles N p in the PIC code is equivalent to NNk, where N 
is the number of velocity grid points and Nk is the number of 
Fourier modes in the continuum code, and the mean squared 
particle weight w 2 (which is called (w 2 ) in Ref.[34il is equiva- 
lent to the continuum value of the mean-square particle distri- 
bution function g 2 . (From the PIC perspective, this is consis- 
tent because the particle weights w are equivalent to Sf/Fo, 
and in w 2 = i w i) /N p , the marker particles have an F 

distribution.) This equivalence between continuum and PIC 
thermal spectra is similar to that found in 2-D hydrodynamics 
between Fourier-Galerkin and point-vortex representations of 
the problem. (The finite-size particles used in most plasma 
PIC codes provides a kind of ultraviolet cutoff that removes 
issues that could arise from point vortices or point particles 
forming tightly-bound pairs.) 

The thermal noise spectrum in PIC codes can be worked 
out using the test-particle superposition principle, assuming 
that shielded test particles can be considered independent and 
random. The equivalence of the PIC and continuum thermal 
spectra indicates that one can likewise consider the random 
phase and amplitude of a Fourier-mode in g at a particular 
velocity (plus the plasma shielding of this mode) to be like 
the random position and weight of a shielded test particle. 

Since the thermal noise level (ip 2 ) = 2^ k (|(p k | 2 ) sca l es 
as g 2 /N, it is important for both PIC and continuum codes to 
either have enough particles or velocity grid points per spa- 
tial grid point so that the noise does not get too large on the 
time scale of the simulation, or to have enough small-scale 
dissipation to prevent the particle weights or g 2 from grow- 
ing too large during the simulation and causing a bottleneck 
problem^ or a numerical diffusion problem!^ (We have con- 
sidered the uniform plasma case in this paper where g 2 is a 
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